1 Creating a model
2 Feature Importance
- 2.1 Mean Decrease Impurity
- 2.2 Permutation Importance
3 Feature Contributions
- 3.1 Plotting feature contributions
4 Joint Feature Contributions
5 PDP and ICE plots
- 5.1 Centered ICE Plots
- 5.2 Two-dimensional PDPs
6 LIME
7 SHAP
8 Resources

In this post I explore different tools and methods that can be used to try and interpret "black box" models.

Since the NFL draft is today, I thought a fun dataset to use for this topic would be the NFL combine data. Using that data we'll first build a model to predict the early career performance for NFL defensive ends (DE). After that, we will use different interpretability methods to better understand the model.

Creating a model¶

Lets get started by importing what we need to set up the data and our modeling pipeline.

In [1]:

%matplotlib inline

import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd

from sklearn.preprocessing import Imputer
from sklearn.model_selection import  cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.metrics import make_scorer
from sklearn.ensemble import RandomForestRegressor

from skll.metrics import spearman

from skopt import BayesSearchCV
from skopt.space import Real, Categorical, Integer

import warnings

In [2]:

# setting up the styling for the plots in this notebook
sns.set(style="white", palette="colorblind", font_scale=1.2, 
        rc={"figure.figsize":(12,9)})
RANDOM_STATE = 420
N_JOBS=8

Here we load in the data that contains each player's combine information and their 3 year Approximate Value, which we will use as measure of early career performance (and going forward will be abbreviated as AV). When building our model we won't be predicting a defensive end's AV, instead we'll be predicting be where his AV ranks (as a percentile) among other defensive ends. Players with high a AV should be closer to 1, while players with low AV should be closer to 0.

NOTE: If you want read up on which combine measurables matter for each NFL position, take a look at Bill Lotter's series of blog posts on the topic.

In [3]:

# read in our data
data_df = pd.read_csv('data/combine_data_since_2000_PROCESSED_2018-04-26.csv')
# onyl get players that have been in the league for 3 years
data_df2 = data_df.loc[data_df.Year <= 2015].copy()
# calculate the player AV percentiles by position
data_df2['AV_pctile'] = (data_df2.groupby('Pos')
                                  .AV
                                  .rank(pct=True,
                                        method='min', 
                                        ascending=True))

Now that we've loaded the data and calculated the AV percentiles, lets get the DE data and create a training set and testing set. We will train and tune our model on the first 8 years (2000-2011) of combine data and then test it on the next 4 years (2012-2015).

In [4]:

# Get the data for the position we want, in this case it's DE
pos_df = data_df2.loc[data_df2.Pos=='DE'].copy().reset_index(drop=True)

# Split the data into train and test sets
train_df = pos_df.loc[pos_df.Year <= 2011]
test_df = pos_df.loc[pos_df.Year.isin([2012, 2013, 2014, 2015])]

We'll be build a random forest model since that seems to be the most common "black box" algorithm people use at work. One thing to to note is that in our modeling Pipeline we will need to include an Imputer because some DEs are missing data as they did not participate in all the combine drills.

In [5]:

# Combine measurables
features = ['Forty',
            'Wt',
            'Ht',
            'Vertical',
            'BenchReps',
            'BroadJump',
            'Cone',
            'Shuttle']
# what we want to predict
target = 'AV_pctile'

X = train_df[features].values
y = train_df[target].values

# the modeling pipeline
pipe = Pipeline([("imputer", Imputer()),
                 ("estimator", RandomForestRegressor(random_state=RANDOM_STATE))])

To tune our model we will use BayesSearchCV from scikit-optimize, which utilizes bayesian optimization to find the best hyperparameters. We'll use Spearman's rank correlation as our scoring metric since we are mainly concerned with the ranking of players when it comes to the draft.

In [6]:

# We use spearman's rank correlation as the scoring metric since
# we are concerned with ranking the players
spearman_scorer = make_scorer(spearman)

# the hyperparamters to search over, including different imputation strategies
rf_param_space = {
    'imputer__strategy': Categorical(['mean', 'median', 'most_frequent']),
    'estimator__max_features': Integer(1, 8),
    'estimator__n_estimators': Integer(50, 500), 
    'estimator__min_samples_split': Integer(2, 200),
}
# create our search object
search = BayesSearchCV(pipe, 
                      rf_param_space, 
                      cv=10,
                      n_jobs=N_JOBS, 
                      verbose=0, 
                      error_score=-9999, 
                      scoring=spearman_scorer, 
                      random_state=RANDOM_STATE,
                      return_train_score=True, 
                      n_iter=75)
# fit the model
# I get some funky warnings, possibly due to the spearman scorer,
# I choose to suppress them
with warnings.catch_warnings():
    warnings.filterwarnings('ignore')
    search.fit(X, y)

In [7]:

# best model parameters
search.best_params_

Out[7]:

{'estimator__max_features': 6,
 'estimator__min_samples_split': 63,
 'estimator__n_estimators': 500,
 'imputer__strategy': 'median'}

In [8]:

# CV score
search.best_score_

Out[8]:

0.38699012444314973

In [9]:

# CV standard deviation
search.cv_results_['std_test_score'][search.best_index_]

Out[9]:

0.13872941403084682

Now that we've tuned our model lets evaluate it on the test set.

In [10]:

# test set data
X_test = test_df[features].values
y_test = test_df[target].values
# predictions
y_pred = search.predict(X_test)
# evaluation
model_test_score = spearman_scorer(search, X_test, y_test)
model_test_score

Out[10]:

0.20929440372118885

So our model's predictions had a Spearman's rank correlation of about 0.209. Is that good or bad? I'm not sure. To get a better sense of how good or how bad of a score that is we can use the actual draft picks as a comparison point.

In [11]:

# create percentiles for nfl draft picks
# Lower numerical picks (e.g. 1, 2, 3) are ranked closer to 1
# Higher numerical picks (e.g. 180, 200, etc) are ranked closer to 0
draft_pick_pctile = test_df.Pick.rank(pct=True,
                                      method='min', 
                                      ascending=False, 
                                      na_option='top')
spearman(y_test, draft_pick_pctile)

Out[11]:

0.6425823905799751

It looks like NFL teams are 3 times better at ranking DEs than our model.

Ok great, now that we have our model lets explore the different tools we can use to better understand it.

Feature Importance¶

Mean Decrease Impurity¶

When using a tree-ensemble like random forest you can find out which features the model found valuable by checking the feature importances. In scikit-learn the feature importances are a reflection of how well a feature reduces some criterion. In our regression example that criterion is mean squared error. This method for calculating feature importance is typically called mean decrease impurity or gini importance.

We can access our model's feature importances with the feature_importances_ attribute.

In [12]:

# get the estimator and imputer from our pipeline, which will be used
# as we try and interpret our model
estimator = search.best_estimator_.named_steps['estimator']
imputer = search.best_estimator_.named_steps['imputer']

estimator.feature_importances_

Out[12]:

array([0.34196584, 0.3615381 , 0.01465109, 0.03990757, 0.05417351,
       0.05790136, 0.07712715, 0.05273537])

Alternatively, we could use eli5's explain_weights_df function, which returns the importances and the feature names we pass it as a pandas DataFrame.

In [13]:

import eli5

# create our dataframe of feature importances
feat_imp_df = eli5.explain_weights_df(estimator, feature_names=features)
feat_imp_df

Out[13]:

	feature	weight	std
0	Wt	0.361538	0.111257
1	Forty	0.341966	0.129695
2	Cone	0.077127	0.106589
3	BroadJump	0.057901	0.088515
4	BenchReps	0.054174	0.081691
5	Shuttle	0.052735	0.079622
6	Vertical	0.039908	0.077301
7	Ht	0.014651	0.037999

So it looks like a player's weight and their forty time are important to model. A thing to note is that explain_weights_df also returns the standard deviations, but they may not be trustworthy as those values assume a normal distribution. Instead of relying on those standard deviations we can access each tree in our ensemble and plot the full distribution of feature importances.

In [14]:

# get the feature importances from each tree and then visualize the
# distributions as boxplots
all_feat_imp_df = pd.DataFrame(data=[tree.feature_importances_ for tree in 
                                     estimator],
                               columns=features)

(sns.boxplot(data=all_feat_imp_df)
        .set(title='Feature Importance Distributions',
             ylabel='Importance'));

Permutation Importance¶

Permutation importances or mean decrease accuracy (MDA) is an alternative to mean decrease impurity that can be applied to any model. The basic idea of permutation importance is to permute the values of each feature and measure how much that permutation negatively impacts the scoring metric (which in our case is the Spearman's rank correlation). This gives us a sense of how our model would perform without that specific feature. All we need to do calculate permutation importance is use PermutationImportance from eli5.

In [15]:

from eli5.sklearn import PermutationImportance

# we need to impute the data first before calculating permutation importance
train_X_imp = imputer.transform(X)
# set up the met-estimator to calculate permutation importance on our training
# data
perm_train = PermutationImportance(estimator, scoring=spearman_scorer,
                                   n_iter=50, random_state=RANDOM_STATE)
# fit and see the permuation importances
perm_train.fit(train_X_imp, y)
eli5.explain_weights_df(perm_train, feature_names=features)

Out[15]:

	feature	weight	std
0	Wt	0.317768	0.038552
1	Forty	0.261409	0.036946
2	Cone	0.033025	0.006485
3	Shuttle	0.023978	0.005698
4	BroadJump	0.017353	0.005326
5	BenchReps	0.016859	0.004783
6	Vertical	0.012590	0.003537
7	Ht	0.005083	0.001254

In [16]:

# plot the distributions
perm_train_feat_imp_df = pd.DataFrame(data=perm_train.results_,
                                      columns=features)
(sns.boxplot(data=perm_train_feat_imp_df)
        .set(title='Permutation Importance Distributions (training data)',
             ylabel='Importance'));

Based on the permutation importances it again looks like Forty and Wt are the two most important features to the model.

Feature Contributions¶

While feature importances can provide us insight into which variables a model finds valuable, they don't tell us how those features impact our model's predictions. One way to actually do that is by using the decision paths in our trees to see how much each feature changes the predictions when going from the parent to the child node. In the end we can break down these feature contributions in a linear manner such that our predictions can be interpreted like so:

$$prediction = bias + feature_{1}contribution + feature_{2}contribution +... feature_{n}contribution$$

Ando Saabas has written a few blog posts that delve into this topic. For now I will try to explain how feature contributions are calculated by looking at an example that uses a shallow tree in our forest.

In [17]:

from IPython.display import Image  
from sklearn.tree import export_graphviz
import graphviz

# Get all trees of depth 2 in the random forest
depths2 = [tree for tree in estimator.estimators_ if tree.tree_.max_depth==2]
# grab the first one
tree = depths2[0]
# now plot the full tree
dot_data = export_graphviz(tree, out_file=None, feature_names=features, 
                           filled=True, rounded=True, special_characters=True)
graph = graphviz.Source(dot_data)  
graph

Out[17]:

Interpretable Machine Learning using Python

Table of Contents

Creating a model¶

Feature Importance¶

Mean Decrease Impurity¶

Permutation Importance¶

Feature Contributions¶